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We describe a Monte Carlo simulation study of the magnetic phase diagram of diluted 
magnetic semiconductors doped with shallow impurities in the low concentration regime. 
We show that because of a wide distribution of interaction strengths, the system exhibits 
strong quantum effects in the magnetically ordered phase. A discrete spin model, found 
to closely approximate the quantum system, shows long relaxation times, and the need 
for specialized cluster algorithms for updating spin configurations. Results for a repre- 
sentative system are presented. 

1. Introduction 

Diluted Magnetic Semiconductors (DMS) doped with a small concentration of (shal- 
low) charged impurities constitute an interesting magnetic system which have a 
number of novel features for study by numerical simulation techniques™! These 
features arise because of the existence of two widely separated length scales - one 
of atomic dimensions (~ 2 A) coming from the d electrons of the magnetic ion, 
and the other from the effective Bohr radius of the shallow impurity, which is typi- 
cally - 20 A. In a certain concentration regime described below, these two lengths 
conspire to produce a very wide distribution of energy scales in the system, which 
poses special challenges for numerical simulations, in particular the problem of long 
equilibration times. 

A prototypical DMS system is a II- VI semiconductor such as CdTe or ZnSe, 
which is in the zinc blende crystal structure, with some of the divalent sites (Cd/Zn) 
substituted by a magnetic ion, most commonly Mn. Manganese has two 4s electrons, 
and so is isovalent with Cd or Zn, but it also has a half filled 3d shell, which by 
Hund's rule, leads to a S — 5/2 ground state. These spins interact via a short range 
exchange coupling in the semiconducting (insulating) host, with nearest and second 
nearest neighbors usually most important. For concentrations x of Mn in, say, 
Cdi-^Mn-jTe, in excess of around 0.2, the Mn spins are connected in a percolated 
network of short range exchange couplings, and the system is found to undergo a 
transition from a paramagnetic state to what is believed to be a spin glass state, at 
a temperature which depends on x, but is of order 10 K. For x well below that, the 
Mn spin system does not percolate, and the magnetic behavior is well described in 
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terms of predominantly isolated Mn spins, and a statistically calculable number of 
small clusters of Mn with typically 2, 3 or 4 Mn ions. 

When shallow dopants are introduced in the system, they are characterized by 
a single electron or hole bound in a hydrogenic Is state, the electron/hole binding 
energy scale being ~ 20 meV. They interact with the Mn d-electron spin via an ex- 
change interaction of the Heisenberg type, which has a typical magnitude ~ 2 meV. 
Thus, it is a reasonable approximation to neglect the modification of the dopant 
electron (or hole") wavefunction by the presence of Mn spins, and the low energy 
Hamiltonian in the limit of low dopant density can be described as a sum of pairwisc 
exchange interactions between the electron spins (sj) and the Mn spins (Sj)B: 

W = ^J(r i ,R J )s i .S j , (1) 

where the sum is over Mn spins at Rj, and electronic spins with wavefunctions 
centered around the donor sites at r^. The exchange interaction energy «7(r$,Rj-) 
between the electron centered at and Mn spin at Rj is proportional to the elec- 
tronic charge density at the Mn site, i.e. proportional to |V>( r i — ) | 2 , i/>( r ) being 
the ground state wavefunction of the dopant electron. (This is very much like the 
contact interaction between nuclear spins and electronic spins, since here the elec- 
tron spin has much greater extent than the spins formed by the Mn 3d electrons.) 
For shallow donors with hydrogenic Is wavefunctions, therefore, 

J(r i ,R j ) = </ e- 2|r *-^ l/aB , (2) 

where Jo is the exchange constant and as the Bohr radius (<~ 20 A). 

For a single dopant, this results in a polarization of the Mn spins at low temper- 
atures in a ferromagnetic configuration with respect to each other in the absence of 
Mn-Mn exchange. For dilute concentration of Mn, so there are only a small fraction 
of near neighbor Mn clusters, we may neglect the direct Mn-Mn interaction, this 
magnetic polaron around each donor electron survives, and in previous work one of 
us has shown! that the relative orientation of the polarons at still lower tempera- 
tures becomes ferromagnetic, as the radius of the polaron increases with lowering of 
T . The system would be expected to show a genuine ferromagnetic transition when 
the magnetic polarons percolate, so that the ferromagnetic order becomes genuinely 
long ranged. 

However, since the percolation fraction for randomly placed spheres in three 
dimensions is ~ 20%, even below the transition, many Mn spins remain essentially 
unattached to the percolated cluster till much lower temperatures. This results 
in a very unusual ferromagnetic phase, in which a substantial part of the spin 
entropy survives down to verv low temperatures. Since analytic results are available 
only in certain limiting caseaJ, we have carried out Monte Carlo simulations of the 

a To simplify the discussion we will consider the electron case where the electronic spin is s = 1/2; 
the hole case with spin-orbit coupling and s = 3/2 is more complex, but not essential for the 
discussion here. Therefore, we will henceforth talk only about the n-doped case. 
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doped DMS system for a typical dopant concentration where the isolated donor 
approximation for the ground state donor wavefunction would be reasonable, and 
a concentration of Mn where the assumption of no direct Mn-Mn interaction would 
be reasonable. 

2. Heisenberg Model 

In most magnetic models involving Heisenberg spins with large quantum numbers 
(S = 5/2 here for the Mn spins which dominate the magnetic response), quantum 
corrections are small and the classical vector spin model has been found to be ade- 
quate. Therefore, we performed Monte Carlo simulations on the model represented 
by Eqs. ([!]) and (g) using classical vector spins. 

On a zinc-blende lattice with lattice constant of 5A, dopants and magnetic ions 
arc distributed randomly on anionic and cationic sites of the lattice respectively. We 
choose cib — 20A, rid — 10 18 cto~ 3 , and x — 0.001. The interaction between electron 
spins and Mn spins is antiferromagnetic. The interaction between electron spins 
has been shown to be unimportant!], while Mn-Mn couplings are not considered 
because of the rare clusters of two Mn spins on the neighboring sites (~ 1%) form a 
magnetically inert singlet at low temperatures. For these concentrations, Mn spins 
are coupled most strongly to only a few electrons, so for each Mn spin we have 
neglected interactions which are smaller than 10~ 5 times its strongest interaction. 
With this cutoff almost all Mn spins are coupled to at least two electron spins. 

We carried out simulations on systems of linear sizes 20, 30, and 40, which have 
256, 864, and 2048 Mn spins and 8, 27, and 64 electron spins, respectivelyu. Since 
we are simulating a very inhomogeneous system with large number of spins, the 
criterion of equilibration is crucial, especially when we try to apply the simulation 
to higher Mn concentrations, when clusters of many nearest neighboring Mn spins 
exist. Therefore, we adopt a widely used scheme in simulating spin glasses to 
determine the equilibration timeEI. Consider two replicas with identical locations of 
particles and identical couplings between them. We start one from a totally random 
spin configuration, and the other from a ferromagnetic configuration, in which Mn 
spins point in the same direction, while electron spins point opposite. We use the 
sample averaged magnitude of magnetization per Mn spin (|M|) as the quantity 
to converge. One expects that increases from zero in the initially random 

replica, whereas it will decrease from the saturation value Mq = 5/2 in the initially 
ferromagnetic replica. After (\M\) of the two replicas agree within error bars, one 
expects both systems have reached equilibrium, which determines the time needed 
before measurements can be taken. 

Figure [l] shows the equilibration of a system of 864 Mn spins at T = 0.00001. 
The replica method turns out to be suitable to determine the equilibration time for 
the Heisenberg model even at a temperature as low as 0.00001. In the simulation, 
we use the standard Metropolis algorithm in which consecutive configurations can 
be different for no more than one single spin. At each Monte Carlo step, we go 
through all the spins successively and flip the spins with equilibrium probabilities. 
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Figure 1: Sample averaged magnetization per spin M/Mq as a function of time 
t for a system of 864 Mn spins at T — 0.00001. The upper curve evolves from 
ferromagnetic initial spin configuration, while the lower one evolves from random 
spin configuration. Classical vector spins are used in the simulation. 
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Figure 2: Normalized magnitude of magnetization per Mn spin (\M\) at different 
temperatures for systems of 256, 864, and 2048 Mn spins in the Heisenberg model. 
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Figure ^] shows the sample averaged magnitude of magnetization per Mn spin 
as a function of temperature for systems of 256 to 2048 Mn spins. (\M\) looks 
very different from that of a conventional ferromagnet. At the lowest temperatures 
simulated for the classical Heisenberg model, (\M\) rises steeply, approaching Mq = 
5/2 in a linear fashion. Magnetization per Mn spin is a good estimate of the 
percentage of Mn spins belonging to the percolating ferromagnetic cluster. The 
critical temperature Tc — 0.014 is obtained by using various cumulants of the 
magnetization, which will be discussed later in Section 3 for the discrete spin model. 




Figure 3: Specific heat C v as a function of temperature for systems of 256, 864, 
and 2048 Mn spins in the Heisenberg model. The solid curve shows the theoretical 
estimate. 

The specific heat per spin is plotted against temperature in Fig. ^. At low 
temperatures, C v approaches unity which is the result of the equipartition theorem. 
This suggests that we cannot treat Mn spins as classical Heisenberg objects. The 
solid curve in Fig. || is a theoretical estimate of C v for the quantum casa-J, which 
is expected to hold at low temperatures (only), and this clearly deviates from the 
simulation results below T c and approaches zero in low temperature limit. It is 
noteworthy that this strong deviation between classical and quantum spins already 
occurs at temperatures where the magnetization is less than 40% of the saturation 
value. Consequently quantum effects are important in this system. 



6 Monte Carlo Simulations of Doped, Diluted Magnetic Semiconductors 



3. Discrete Spin Model 



3.1. Discrete 12- state model 

In order to simulate the effects of the discreteness of the quantum spins in doped 
diluted magnetic semiconductors without having to perform a fully quantum simu- 
lation, we have modified Mn spins from a continuous vector spin to various discrete 
spin versions. Here, we discuss the results of a discrete 12-state model, in which each 
spin can take one of the 12 (110) directions of a cubic lattice with a Hamiltonian 
given by Eqs. ([!]) and (g). 

Figure |(a) and (b) compare the evolution of magnetization per Mn spin in a 
system of 256 Mn spins at T — 0.04 and T = 0.00001 using standard single spin 
flip dynamics. At the higher temperature, of the two replicas converge to 

the same value, implying the equilibrium is reached in both replicas. However, 
at temperature three orders below, the magnetization values hardly change after 
initial 100 steps and show no indication of equilibration. This is easily understood as 
follows: each single spin flip is associated with a finite energy change. Therefore, the 
discretization of the energy levels will inevitably lead to significant slowing down of 
the simulation at low temperatures. The broad distribution of interactions worsen 
the equilibration of discrete spins by forming domains around each electron spins 
where the exchanges are strong. These domains are unlikely to flip in reasonable 
times, because the exchanges at the domain surfaces are much weaker compared 
to that in the bulk. A new algorithm, which can effectively flip each domains at a 
single step, is thus required to achieve equilibration. 

3.2. Cluster algorithm 

To develox) a cluster algorithm for various spin models, we follow Wolff's cluster 
algorithm!! We first define Ising-like spin-flip operation s — > s' as the reflection 
with respect to the plane orthogonal to a unit vector n, 



which satisfies 
and 



1Z(n)s = s — 2(s ■ n)n, 
n{n) 2 = 1, 

[ft(n)si] • [K{h)s 2 ] = si • s 2 . 



(3) 

(4) 
(5) 



In the discrete spin model, only finite number of unit vectors n exist that transform 
an arbitrary discrete spin into other discrete spins. For instance, unit vectors along 
(110) and (100) directions are a set of appropriate reflections for discrete spins along 
(110) directions. 

Due to the non-nearest-neighboring nature of spin exchanges in doped magnetic 
semiconductors, the cluster algorithm needs extra consideration whether the whole 
cluster is flipped or not. The revised cluster algorithm now consists of the following 
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Figure 4: Magnetization per spin M as a function of time t for system with 
256 Mn spins in replicas starting from ferromagnetic configuration and random 
configuration in the discrete 12-spin model (a) at T = 0.04 with single spin flips, 
(b) at T = 0.00001 with single spin flips, and (c) at T = 0.00001 with the cluster 
updating algorithm described in the text. 
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sequence of operations, (a) Choose a electron spin So at ro and an appropriate 
reflection n randomly. So is the seed of the cluster C to be built, (b) Visit all 
manganese spins and activate the bond between the electron spin and a manganese 
spin S at R with probability 

Vi(S) = l-cxp{min[0,2J(ro,R)(s -n)(S-n)/T]}. (6) 

If the bond is activated, add S to the cluster C. (c) Visit all electron spins, except 
the seed so, and add an electron spin s to the cluster C with probability 

V n (s) = 1 — exp{min[0, ^ 2J(r, R)(s ■ h)(S ■ h)/T]}. (7) 
sec 

(d) Flip the whole cluster C with probability 

Vin(C) = exp{min[0, 2J(r,R)(s • n)(S • h)/T]}. (8) 

sec-{s },sec 

Detailed balance is fulfilled since the transition probabilities W between two con- 
figurations {s, S} and {s', S'} that differ by a flip 7Z(n) on a cluster C built around 
electron spin Sq obey 

W({s,S}^{s',S'}) _ Pz/z(C)n se c[l-^(s)]n sSC [l-^(S)] 



W({s', S'} - {s, S}) V in {C>) n., gc ,[l - Vh( S ')} Us^a I 1 - W)] 



exp 



exp 




sec, sec s=s ,sec 



£ J(r,R)(s-n)(S-n) 



-i^J(r,R)( S '.S'- S -S)^ (9) 



;.S 



It is worth pointing out that all probabilities for activating bonds within C are the 
same before or after the cluster flipping because of Eqs. (^) and (Q). Due to the 
extremely broad distribution of spin exchanges, we add the following step to reduce 
the time in which all the manganese spins far away from any electron spin can reach 
thermal equilibrium, (e) Visit all manganese spins and flip each individual spin S 
to an arbitrary discrete spin S' with probability 

V(S) = exp{min[0, ^ ^( r . R ) s ■ ( s - S')/T\} (10) 

s 

Ergodicity of the above processes is guaranteed by the fact that there is always 
nonzero probability that C consists of only one spin, and that there is always one 
reflection or two consecutive reflections connecting any two discrete spins. This 
allows exploration of all configurations with finite probabilities required for ergod- 
icity. 

The efficiency of the cluster algrithm can be tested at very low temperatures. 
Figure ||(c) shows the average evolution curve of the systems of 256 Mn spins at 
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T = 0.00001. Cluster algorithm is applied in the random replica. In 100,000 steps, 
(\M\) in the random replica rises to almost 97% of the value in the ferromagnetic 
replica. The small discrepancy largely comes from the worst 10% of the samples in 
which a longer time is required to remove all the domain walls. 

The development of ordered state from initially random configuration particu- 
larly in systems of 256 Mn spins using cluster algorithm confirms that the replica 
with such ferromagnetic initial configuration can reach equilibrium by single spin- 
flip dynamics. Obviously, the reason is that, for the simple Hamiltonian in Eq. (jl]), 
we know that the ground state is a ferromagnetic state at zero temperature, since 
there is no frustration. Therefore, even if the random replica is not in equilib- 
rium, we can measure the equilibrium quantities in the replica starting from the 
appropriate initial configuration after enough many steps using only single spin-flip 
algorithm. 

Although the cluster algorithm with the replica method determining the equi- 
libration time may appear to be a luxury in the dilute concentration limit, it is 
absolutely necessary when sizable clusters with Mn-Mn interactions exist at higher 
Mn concentrations when the ground state configuration is unknown. 

3.3. Results of the discrete 12-state model 

In this subsection, we show results based on the assumption that the replica starting 
from ferromagnetic initial configuration reaches equilibrium after long enough time, 
using single spin-flip dynamics. 

A dimensionless quantitycl used to characterize the ferromagnetic order is given 



The coefficients in Eq. ( |l l| ) are chosen so that Gl approaches zero in the thermody- 
namic limit above T c and unity below T c . From finite size scaling theory, we expect 
that Gl is size independent at the critical temperature. The three Gl curves in 
Fig. H seem to cross at T c ~ 0.03. Below T c Gl increases as system size increases, 
while Gl decreases as system size increases above T c , as expected. 

The magnetization per Mn spin is plotted versus temperature for the three 
different system sizes in Fig ^. Compared with the results in the Heisenberg model, 
the curves are roughly scaled by a factor of 2 on the temperature axis, as is T c . The 
specific heat per spin shown in Fig. ^ starts dropping at T ~ 0.001. The simulation 
results in the discrete 12-state model agree xmite well with the theoretical estimate 



by 




(11) 
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Figure 6: Normalized magnitude of magnetization per Mn spin versus temperature 
for different system sizes in the discrete 12-state model. The solid line fits (\M\) of 
the system of 864 Mn spins in the Heisenberg model. 
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4. Conclusion 

In conclusion, we have presented results of a numerical simulation appropriate for 
diluted magnetic semiconductors with a low concentration of magnetic ions (x) and 
a low concentration of shallow dopants {rid). The restriction on rid is provided by the 
requirement that the system remain insulating, i.e. below the Mott concentration 
(or explicitly rid • a% < 0.016). The restriction on x is required for the Mn-dopant 
interaction to be the dominant interaction. 

In this regime, the system is found to order ferromagnetically at low tempera- 
tures, but the ferromagnetic state is rather unusual, in that the magnetization does 
not approach the saturation value till well below the ordering temperature. This 
arises because of the effective couplings that characterize the system are distributed 
on a logarithmic scale owing to the exponential variation of the interactions with 
distance. Because of this wide distribution, the system cannot be adequately rep- 
resented by a classical Heisenberg model at the low temperatures of interest. A 
discrete spin model is found to remove most of the inadequacies of the continuous 
spin model, but requires the application of specialized cluster algorithms to help 
alleviate the extremely long relaxation times that are inherent in a discrete spin 
model. 

While there may be other methods to get around the problems of equilibra- 
tion in the ferromagnetic regime, the cluster algorithms described above, which are 
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found to work well (at least for small sizes) in the concentration regimes studied, 
should in principle be useful in other regimes, e.g. higher Mn concentration, where 
Mn-Mn interactions cannot be neglected. The higher concentration regimes will 
involve frustration effects as well, and having a reliable numerical method to ensure 
equilibration is prerequisite for any study of such effects in the diluted magnetic 
semiconductors. We are currently testing if big enough sizes can be brought to 
equilibrium with the above scheme to enable finite size scaling for accurate location 
of T c , and also to ensure that finite size corrections for the susceptibility, specific 
heat and other thermodynamic quantities in the ferromagnetic phase remain small. 
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